clear;
load('mdc_doping.mat');

global mdc

sidepr=[];
mpp=[];%major peak positions
spp=[];%side peak positions
er=[];%error
figure('position',[300,300,200,500]); hold on; box on;
for ii=1:length(MDC);

mdc=MDC{ii}*0.5e5;
kk=K{ii};

mdc=mdc-min(mdc);
mdc=mdc(abs(kk)<0.788);
mdc=(mdc+flipud(mdc))/2;
kk=kk(abs(kk)<0.788);
mdc=mdc./max(mdc)*0.5;


if isempty(mpp)
    x0=[0.16,0.2,0.12,0.05,0.5,0.1];
else
    x0=x;
end


[x,fval,exit,out,lamb,grad,hessian]=fmincon(@FindMDCPeak,x0,[],[],[],[],[0,0,0.01,0.001,0.47,0.01],[1,1,1,1,0.7,0.5]);



A1=x(1);%amplitude of peak 1
m1=x(2);%mean of peak 1
g1=x(3);%gamma of peak 1
A2=x(4);
m2=x(5);
g2=x(6);

%estimate error
dA1=sqrt(FindMDCPeak(x)/hessian(1,1));
dA2=sqrt(FindMDCPeak(x)/hessian(4,4));
dr=dA2/A1+A2*dA1/A1^2;

k=linspace(1,length(mdc),length(mdc));
k=k-mean(k);
k=k/max(k)*0.788;

L1=A1*(g1^2./((k-m1).^2+g1^2))/(pi*g1);
L2=A2*(g2^2./((k-m2).^2+g2^2))/(pi*g2);
L=L1+L2;
L=(L+fliplr(L));

offset=3.3;%y axis offset
inteval=0.53;%y axis offset
plot(k,0*k-inteval*ii+offset,'k');
plot(kk,mdc-inteval*ii+offset,'r');
plot(k,L-inteval*ii+offset,'b');
plot(k,L1-inteval*ii+offset,'g');
plot(k,L2-inteval*ii+offset,'c');
plot(k,fliplr(L1)-inteval*ii+offset,'m');
plot(k,fliplr(L2)-inteval*ii+offset,'y');
% 

sidepr=[sidepr,max(L2)/max(L1)];
mpp=[mpp,m1];
spp=[spp,m2];
er=[er,dr];
end
xlim([-0.788,0.788]);
ylim([0,3.4]);
xlabel('k_{||} (�)');
ylabel('intensity (arb. units)');
set(gca,'PlotBoxAspectRatio',[1 3 1],'FontSize',12);

figure('position',[100,100,400,300]);
errorbar(doping,sidepr,er*2,'o','MarkerSize',10);
xlim([0,45]);
set(gca,'xtick',0:10:40)
xlabel('doping(%)');
ylabel('folding peak relative intensity')
ylim([0,0.4])
set(gca,'PlotBoxAspectRatio',[3 1 1],'FontSize',12);
